Author

Emi Carlevaro

Published

January 22, 2026

tHIS IS THE SIMULATION i HAD IN THE stc PROJECT. BUT TI GIVES poor identification of parameters, with large sets. Even when using const',HC3orNeweyWest` estimators.

Simulation data comes from C:\Users\Emi\OneDrive\UWA PhD\MonPolicy StockMkt\Results\Simulations\NoCommShockZ\CalCrisis-HARRVBig-NoZ-1016-Tres\simVARs.RDS[[1]]

see details in that folder.

INIT

Code 
#required packages to make a package
pacman::p_load(devtools, roxygen2, testthat)
#dependecies
#pacman::p_load(tidyverse, data.table, DT, plotly, gridExtra, latex2exp, here, logr, future.apply, furrr)
#pacman::p_load(dplyr,parallel,sandwich,here,future,logr,profvis,lobstr,parallel)
#https://ourcodingclub.github.io/tutorials/writing-r-package/
#create_package
#https://combine-australia.github.io/r-pkg-dev/functions.html
#load_all()
#run these two lines and then library(GridMaker) in normal rstudio to install package
# load_all loads a package. It roughly simulates what happens when a package is installed and loaded with library().
devtools::load_all()

devtools::install()
Code 
library(vars)
library(tidyverse)
library(plotly)
library(iStabi)
Code 
sam <- iStabi::data_sim_biVAR

0. THE MODEL

The simulated model is a bivariate SVAR(2) as follows \[ \begin{equation*} \begin{bmatrix} 1 & -\beta \\ -\alpha & 1 \end{bmatrix} \begin{bmatrix} i_t \\ s_t \end{bmatrix} = \begin{bmatrix} 0.3 & 0.3 \\ -0.2 & 0.6 \end{bmatrix} \begin{bmatrix} i_{t-1} \\ s_{t-1} \end{bmatrix} + \begin{bmatrix} -0.2 & 0.1 \\ 0.1 & 0.1 \end{bmatrix} \begin{bmatrix} i_{t-2} \\ s_{t-2} \end{bmatrix} + \begin{bmatrix} \varepsilon_t \\ \eta_t \end{bmatrix} \end{equation*} \] where the impact coefficients are set similar to our findings for the before-ELB sample: \(\alpha=-5,\;\beta \times 100 = 1\) . The error terms \(\left\{ \varepsilon_{t}, \eta_{t}, \right\} _{t=1}^{T}\) are independently normally distributed with mean zero and time-varying variances \(\sigma^{2}_{\varepsilon; t}\), \(\sigma^{2}_{\eta; t}\).

1. ESTIMATION VAR

Code 
varModel <- VAR(cbind('i' = sam$i, 's' = sam$s), p = 2, type = "const")
Res = residuals(varModel)
Code 
# Use for computing the moment function
Omega <- cbind('ii' = Res[, 'i']*Res[, 'i'], 'is' = Res[, 'i']*Res[, 's'],
               'ss' = Res[, 's']*Res[, 's'])

Y <- cbind( 'omega_ii' = Omega[, 'ii'], 'omega_ss' = Omega[,'ss'], 'omega_is' = Omega[,'is'] )

# Compute the moment function at each t. 
# INPUT: a named vector with values for the structural parameters
# OUTPUT: a TxG matrix where each row corresponds to the value of the moment function on that obseration (date)
# an each column the value fot ehfat moment equation
GET_MOM_VALUES <- function(t0, ...) {
  list2env(as.list(t0),envir = environment())
  
  #Y <- cbind( 'omega_ii' = Res[,'i']^2, 'omega_ss' = Res[,'s']^2, 'omega_is' = Res[,'i'] * Res[,'s'] )
  
  #Y %*% c(-theta0['alpha'],
  #        -theta0['beta'],
  #        1 + theta0['alpha'] * theta0['beta'])
  
  #fT 
  as.matrix(-alpha*Omega[, 'ii'] - beta*Omega[, 'ss'] + (1+alpha*beta) * (Omega[, 'is']))

  
}
Code 
GET_B <- function(theta0) {
  c(-theta0['alpha'],
    -theta0['beta'],
    1 + theta0['alpha'] * theta0['beta'])
}

3. BUILD SETS

Critical values

Critical values for \(qLL-\tilde{S}\) for 1 moment condition from table VII Suplementary material, p3 Magnusson Mavroeidis (2014).

THe S-test follows a Chi-Square distribution with 1 degrees of freedom (equatl to the number of moment conditions).

Grid has 520 rows.

The number of fixed identified parameters is 0.

The number of estimated parameters under the null is 0.

Code 
NSIPARAMS = 0 # not needed in the new version of the package
PCONFS = c(0.80, 0.85)
get_CVs <- function() {
  
  CVs <- list('qLLStab'=NA, 'S'=NA, 'genS_qLL'=NA)

  CVs$qLLStab <- filter(iStabi::data_CVs$qLLStab, DG_F == KZ & P_CONF %in% PCONFS)
  CVs$genS_qLL <- filter(iStabi::data_CVs$genS_qLL, 
                         MOMENT_CONDITIONS == KZ & STRONGLY_IDENTIFIED_PARAMETERS == NSIPARAMS & P_CONF %in% PCONFS)
  CVs$S <- tibble('P_CONF' = PCONFS,
                 'CUTOFF_TOP' = map_dbl(PCONFS, ~qchisq(.x, KZMKX)))
  
  CVs
}

CVs <- get_CVs()

Critical values are:

Code 
CVs$S
# A tibble: 2 × 2
  P_CONF CUTOFF_TOP
   <dbl>      <dbl>
1   0.8        1.64
2   0.85       2.07
Code 
CVs$qLLStab
# A tibble: 2 × 3
   DG_F P_CONF CUTOFF_TOP
  <dbl>  <dbl>      <dbl>
1     1   0.8        5.92
2     1   0.85       6.46
Code 
CVs$genS_qLL
# A tibble: 2 × 4
  MOMENT_CONDITIONS STRONGLY_IDENTIFIED_PARAMETERS P_CONF CUTOFF_TOP
              <dbl>                          <dbl>  <dbl>      <dbl>
1                 1                              0   0.8        7.14
2                 1                              0   0.85       7.74

Data for plotting

Code 
build_set <- function(setName) {
  thisCV <- CVs[[setName]]
  statMax = filter(thisCV, P_CONF == PCONFS[2])$CUTOFF_TOP
  thisPconf = thisCV$P_CONF
  thisCV = thisCV$CUTOFF_TOP

  theseCols <- c(PARAMS_CFG$name, setName)
  aSet = Grid[get(setName) <= statMax, ..theseCols][, 
                                                    sigAt := thisPconf[ 1+findInterval(get(setName), thisCV,left.open=TRUE)]]
  aSet
  
}
Code 
Sets <- list('qLLStab' = build_set('qLLStab'),
             'S' = build_set('S'),
             'genS_qLL' = build_set('genS_qLL'))
saveRDS(Sets, 'Sets_Unidentified_model.rds')

4. PLOTTING

Grid for plotting

Force the starting values to be a multiple of the step size. This avoids problems and is clear what the true parameter space is.

Plots

True parameter vector $ =(= -5, = 0.01)$

genS-qLL plot

Code 
# INTERACTIVE PLOTs ----
pSbase <- plot_ly(data = Sets$genS_qLL[sigAt<=0.80,], 
                 x=~alpha, 
                 y=~beta,
                 color='grey')%>%
  layout(
    title = "genS_qLL (80% confidence)",
    xaxis = list(
      title = plotly::TeX("\\alpha"),
      range = c(-20, 20),
      fixedrange = TRUE
    ),
    yaxis = list(
      title = plotly::TeX("\\beta"),
      tickvals = seq(from = -0.03, to = 0.03, by = 0.005),
      range = c(-0.02, 0.02),
      fixedrange = TRUE
    )
  ) %>%
  add_trace(
    x = c(-5),
    y = c(0.01),
    type = "scatter",
    mode = "markers",
    marker = list(symbol = "square", size = 12, color = "red"),
    inherit = FALSE,
    showlegend = FALSE
  ) %>%
  config(mathjax = "cdn")

pSbase %>%
  add_trace(type='scatter', mode='markers', marker = list(color='blue'))
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
Warning in RColorBrewer::brewer.pal(N, "Set2"): minimal value for n is 3, returning requested palette with 3 different levels

Static plot: qLL-Stab

Code 
pBase <- ggplot(data=Sets$qLLStab[sigAt<=0.80,], 
                mapping=aes(x=alpha, y=beta)) +
  scale_x_continuous(name=TeX(r"(\\alpha)"), limits=c(-20, 20), 
                     breaks=seq(from=-20, to=20, by=2)) +
  scale_y_continuous(name=TeX(r'($\\beta$)'), limits=c(-0.03, 0.03), 
                     breaks=seq(from=-0.03, to=0.03, by=0.01)) + 
  theme_bw() + theme(panel.grid=element_line(colour='#999999', linetype='14'), 
                     panel.grid.minor=element_line(colour='white'), 
                     text=element_text(size=10), 
                     axis.title.x = element_text(colour='red'),
                     axis.title.y = element_text(colour='blue', size=rel(0.9), hjust=0.9),
                     axis.text.x = element_text(size=rel(1.25), angle=90, vjust=0.5), 
                     axis.text.y = element_text(size=rel(1.25))) +
  geom_vline(xintercept=0, colour = 'black', linetype = 'solid') +
  geom_hline(yintercept=0, colour = 'black', linetype = 'solid') +   
  # Highlightinh the TRUE value
  geom_point(aes(x=-5, y=0.01), size=3, colour='black', fill='yellow')

# Scatterplot
pGplot <- pBase + geom_point(aes(colour=qLLStab))
ggsave('Unidentified_model_qLLStab_set.png', plot=pGplot, width=6, height=5)
Warning in geom_point(aes(x = -5, y = 0.01), size = 3, colour = "black", : All aesthetics have length 1, but the data has 155 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.
Code 
pGplot
Warning in geom_point(aes(x = -5, y = 0.01), size = 3, colour = "black", : All aesthetics have length 1, but the data has 155 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing
  a single row.